library(TCGAbiolinks)
library(tidyverse)
library(tidySummarizedExperiment)
library(survival)
library(survminer)

read_delim('00_util_scripts/ref/hg38_gtf-ensembl-symbol.csv')

slice_expr_quart <- function(df, gene){
  df |>
    mutate(
      goi = {{ gene }},
      vital_status = ifelse(vital_status == 'Alive', 0, 1),
      expr = case_when(goi <= quantile(goi, c(.25)) ~ 'lowest 25%',
                       goi >= quantile(goi, c(.75)) ~ 'highest 25%',
                       .default = 'NA')) |>
    filter(expr != 'NA')
}

# CTSP-DLBCL1 ----------
query <- GDCquery(project = 'CTSP-DLBCL1',
                  data.category = 'Transcriptome Profiling',
                  data.type = 'Gene Expression Quantification')

# GDCdownload(query, files.per.chunk = 1)

## default summarizedExpr object integration always fail for this proj
CTSPdata <- GDCprepare(query, summarizedExperiment = FALSE) |>
  select(gene_name, contains('fpkm_uq')) |>
  as_tibble()

CTSPdata |>
  write_rds('mission/DLBCL-Btk/data/CTSP-DLBCL.rds')

CTSPdata <- read_rds('mission/DLBCL-Btk/data/CTSP-DLBCL.rds')

# manual download clinical data
query <- GDCquery(project = 'CTSP-DLBCL1',
                  data.category = 'Clinical')

# GDCdownload(query, files.per.chunk = 1)

# write a func to parse the fucking too-free json
parse_ctsp <- function(df){
  df |>
    mutate(vital_status = discard(vital_status, is.na),
           gender = discard(gender, is.na),
           days_to_follow_up = max(days_to_follow_up, na.rm = TRUE),
           alt_id = submitter_id[[3]] |> str_remove('-.+')) |>
    slice_head(n = 1)
}

ctsp_meta <- list.files(path = 'GDCdata/CTSP-DLBCL1/', pattern = '.json', recursive = TRUE, full.names = TRUE) |>
  map(jsonlite::fromJSON) |>
  map(pluck, "ClinicalData", "GDC Clinical Data Entities") |>
  map(select, c('submitter_id', 'vital_status', 'days_to_follow_up', 'gender')) |>
  map(parse_ctsp) |>
  list_rbind()
  
gene_ctsp <- CTSPdata |>
  filter(gene_name %in% c('MARCKS')) |>
  pivot_longer(where(is.numeric), names_to = 'sample', values_to = 'expr') |>
  pivot_wider(names_from = gene_name, values_from = expr) |>
  mutate(sample = str_remove_all(sample, '.+stranded_|-sample|-TTP.+'))

ortho_meta <- gene_ctsp |>
  filter(str_detect(sample, 'CTSP')) |>
  left_join(ctsp_meta, join_by(sample == submitter_id))

ctsp_final <- gene_ctsp |>
  filter(!str_detect(sample, 'CTSP')) |>
  left_join(ctsp_meta, join_by(sample == alt_id)) |>
  filter(!is.na(vital_status)) |>
  bind_rows(ortho_meta) |>
  mutate(sample = case_when(str_detect(sample, 'DLBCL') ~ submitter_id,
                            .default = sample)) |>
  select(-c(submitter_id, alt_id))

ctsp_final |>
  write_csv('mission/DLBCL-Btk/results/CTSP-DLBCL-survival.csv')

ctsp_final |>
  slice_expr_quart(MARCKS) |>
  survfit(Surv(days_to_follow_up, vital_status) ~ expr, data = _) |>
  ggsurvplot(pval = TRUE,
             risk.table = 'nrisk_cumcensor',
             title = "CTSP-DLBCL patients OS with MARCKS expression")

# TCGA-DLBC -------------
query <- GDCquery(project = 'TCGA-DLBC',
                  data.category = 'Transcriptome Profiling',
                  data.type = 'Gene Expression Quantification')

# api download can be interrupted frequently and very slow
# recommand set files.per.chunk = 1 or gdc-client from gdc website
# GDCdownload(query, files.per.chunk = 1)

# GDCprepare() can be useful to gather
DLBCdata <- GDCprepare(query)

write_rds(DLBCdata, 'DLBCL-Btk/data/TCGA-DLBC.rds')

DLBCdata <- read_rds('DLBCL-Btk/data/TCGA-DLBC.rds')

gene_expr <- DLBCdata |>
  assay('fpkm_uq_unstrand') |>
  as_tibble(rownames = 'ENSEMBL') |>
  mutate(ENSEMBL = str_remove(ENSEMBL, '\\..+')) |>
  right_join(gene_list) |>
  select(-ENSEMBL) |>
  pivot_longer(where(is.numeric), names_to = 'sample', values_to = 'fpkm_uq') |>
  pivot_wider(names_from = 'SYMBOL', values_from = 'fpkm_uq')

dlbc_meta <- DLBCdata@colData |>
  as_tibble() |>
  select(barcode, patient, vital_status, days_to_death, days_to_last_follow_up, gender, age_at_index) |>
  mutate(days = case_when(!is.na(days_to_death) ~ days_to_death,
                          .default = days_to_last_follow_up),
         .keep = 'unused') |>
  right_join(gene_expr, join_by(barcode == sample)) |>
  select(-barcode)

dlbc_meta |>
  write_csv('DLBCL-Btk/results/TCGA-DLBC-survival.csv')

dlbc_meta |>
  slice_expr_quart(MARCKS) |>
  survfit(Surv(days, vital_status) ~ expr, data = _) |>
  ggsurvplot(pval = TRUE,
             risk.table = 'nrisk_cumcensor',
             title = "TCGA-DLBC patients OS with MARCKS expression")

## unrealistic to directly merge two datasets
dlbc_meta |> ggplot(aes('',MARCKS)) + geom_violin()
ctsp_final |> ggplot(aes('',MARCKS)) + geom_violin()

# NCICCR-DLBCL -----------
query <- GDCquery(project = 'NCICCR-DLBCL',
                  data.category = 'Transcriptome Profiling',
                  data.type = 'Gene Expression Quantification')

# api download can be interrupted frequently and very slow
# recommand set files.per.chunk = 1 or gdc-client from gdc website
# GDCdownload(query, files.per.chunk = 1)

# GDCprepare() can be useful to gather
NCICCRdata <- GDCprepare(query)

assays(NCICCRdata) <- assays(NCICCRdata)[c("unstranded","fpkm_uq_unstrand")]

write_rds(NCICCRdata, 'DLBCL-Btk/data/NCICCR-DLBCL.rds')

NCICCRdata <- read_rds('DLBCL-Btk/data/NCICCR-DLBCL.rds')

nci_expr <- NCICCRdata |>
  assay('fpkm_uq_unstrand') |>
  as_tibble(rownames = 'ENSEMBL') |>
  mutate(ENSEMBL = str_remove(ENSEMBL, '\\..+')) |>
  right_join(gene_list) |>
  select(-ENSEMBL) |>
  pivot_longer(where(is.numeric), names_to = 'sample', values_to = 'fpkm_uq') |>
  pivot_wider(names_from = 'SYMBOL', values_from = 'fpkm_uq')

nci_meta <- NCICCRdata@colData |>
  as_tibble() |>
  select(sample, vital_status, days_to_last_follow_up, gender, age_at_diagnosis) |>
  filter(!is.na(days_to_last_follow_up)) |>
  left_join(nci_expr)

nci_meta |>
  write_csv('DLBCL-Btk/results/NCICCR-DLBCL-survival.csv')

nci_meta |>
  slice_expr_quart(MARCKS) |>
  survfit(Surv(days_to_last_follow_up, vital_status) ~ expr, data = _) |>
  ggsurvplot(pval = TRUE,
             risk.table = 'nrisk_cumcensor',
             title = "NCICCR-DLBCL patients OS with FCGR2B expression")

# sum up ------
tcga_meta <- read_csv('DLBCL-Btk/results/TCGA-DLBC-survival.csv') |>
  slice_expr_quart(MARCKS)

ctsp_meta <- read_csv('DLBCL-Btk/results/CTSP-DLBCL-survival.csv') |>
  rename(days = days_to_follow_up) |>
  slice_expr_quart(MARCKS)

nci_meta |>
  rename(days = days_to_last_follow_up) |>
  slice_expr_quart(MARCKS) |>
  bind_rows(tcga_meta, ctsp_meta) |>
  mutate(vital_status = if_else(vital_status == 'Dead', 1, 0)) |>
  survfit(Surv(days, vital_status) ~ expr, data = _) |>
  ggsurvplot(pval = TRUE,
             risk.table = 'nrisk_cumcensor',
             legend.labs = c('Highest 25%','Lowest 25%'),
             title = "TCGA-DLBCL patients OS with MARCKS expression") +
  labs(x = 'days')
